/*-----------------------------------------------------------------------------------
Name: Sangyoon Park
Date: February 4 2023
This do file does : 
	Controls False Discovery Rate to adjust significance levels
-----------------------------------------------------------------------------------*/

clear all
set matsize 10000
set maxvar 10000
set more off
capture log close


#delim ;

/*-----------------------------------------------------------------------------------
			 Set up
-----------------------------------------------------------------------------------*/

		
/*-----------------------------------------------------------------------------------
			 		DATA = GAME
-----------------------------------------------------------------------------------*/		

	import excel "${raw_path}/DATA_GAME_ALL.xlsx", sheet("Sheet1") firstrow ;
	tempfile game;
	save `game';
	clear;		
	
/*-----------------------------------------------------------------------------------
			 		DATA = EXPORTER SURVEY
-----------------------------------------------------------------------------------*/		

	use "${data_path}/INTERMEDIARY_LONG.dta", clear;
	rename age firmage;
		egen group_id = group(strata st_tgroup);
		gen tr_farm = (st_tgroup == 2 | st_tgroup == 6);
		gen tr_int = (st_tgroup == 3 | st_tgroup == 7);
		gen tr_joint = (st_tgroup == 4 | st_tgroup == 8);
		gen certelig = (st_tgroup == 5 | st_tgroup == 6 | st_tgroup == 7 | st_tgroup == 8);
		gen tr_fXcert = tr_farm * certelig;
		gen tr_iXcert = tr_int * certelig;
		gen tr_jXcert = tr_joint * certelig;	
	tempfile intdata;
	save `intdata';
	clear;	

/*-----------------------------------------------------------------------------------
			 		DATA = FARMER SURVEY
-----------------------------------------------------------------------------------*/		


*** read data;

	use "${data_path}/FARMER_LONG_SETUP", replace;
	
*** Drop control communes;

	drop if control_commune == 1;	
	
*** Generate dummy for treatment groups;	
	
	tab st_tgroup, gen(dum_tgroup);
	
*** Cluster for standard errors;

	egen group_id = group(strata st_tgroup);

*** Generate treatment levels and interaction terms;

	gen tr_farm = (st_tgroup == 2 | st_tgroup == 6);
	gen tr_int = (st_tgroup == 3 | st_tgroup == 7);
	gen tr_joint = (st_tgroup == 4 | st_tgroup == 8);
	gen certelig = (st_tgroup == 5 | st_tgroup == 6 | st_tgroup == 7 | st_tgroup == 8);
	gen tr_fXcert = tr_farm * certelig;
	gen tr_iXcert = tr_int * certelig;
	gen tr_jXcert = tr_joint * certelig;

*** Drop large farms (farm size larger than five hectares);

	drop if ha_whitedragon >= 5;	
	
*** Define local variables;
	
	local dummy_coef dum_tgroup2 dum_tgroup3 dum_tgroup4 dum_tgroup5 dum_tgroup6 dum_tgroup7 dum_tgroup8;
	local treat_coef "tr_farm tr_int tr_joint certelig tr_fXcert tr_iXcert tr_jXcert";

	local baseline "age female edu_second exp_dragon ha_whitedragon tree_whitedragon anycert_base loan_any_base save_bank_base mean_trust_base mean_business_base mean_confidence_base raven_total present_bias_base gap_comply_base hours_base totalvolume_base price_base cost_fert_base cost_pest_base cost_facequip_base cost_labor_base cost_utility_base cost_total_base attrit_fu1 attrit_fu2";
	local balance age female edu_second exp_dragon ha_whitedragon white_age_less1 white_age_less3 white_age_less10 white_age_less20 white_age_more20 anycert_base loan_any_base save_bank_base mean_trust_base mean_business_base mean_confidence_base raven_total present_bias_base gap_comply_base hours_base totalvolume_base price_base log_cost_fert_base log_cost_pest_base log_cost_facequip_base log_cost_labor_base log_cost_utility_base;
	local intbalance int_age int_size_facility int_volume_china_fu0 int_volume_asia_fu0 int_volume_eu_fu0 int_volume_dom_fu0 int_dum_firmtype1 int_dum_firmtype2 int_dum_firmtype3;

	
	tempfile maindata;
	save `maindata';
	clear;	

		
/*-----------------------------------------------------------------------------------
			 		To combine with FDR q-values
-----------------------------------------------------------------------------------*/	

/*-----------------------------------------------------------------------------------


			 		Tables for Main Text

					
-----------------------------------------------------------------------------------*/	


/*-----------------------------------------------------------------------------------
			 		OUTCOME = GAP Audit by Management Area (Table 2)
-----------------------------------------------------------------------------------*/	


	use `maindata';
	
	local outcome "score pestman equipman hygman landman fertman";
	local count_outcome = `:word count `outcome'';
	mat X = J(7*`count_outcome',1,0);
	local i = 1;
	
	foreach y of local outcome {;
	reg std_`y' tr_farm tr_int tr_joint certelig tr_fXcert tr_iXcert tr_jXcert `balance' `intbalance' i.strata i.round, robust cluster(group_id);

	foreach z of local treat_coef {;
	local t_score = (_b[`z']/_se[`z']);
	local p_score = (2*ttail(e(df_r),abs(`t_score')));
	mat X[`i',1] = `p_score';
	local i = `i' + 1;
	};
	};
	
	clear;
	
*** Anderson's Code for FDR q-value;	
	
	svmat X;
	rename X1 pval;
	
* Collect the total number of p-values tested;

quietly sum pval;
local totalpvals = r(N);

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank;

quietly gen int original_sorting_order = _n;
quietly sort pval;
quietly gen int rank = _n if pval~=.;

* Set the initial counter to 1 ;

local qval = 1;

* Generate the variable that will contain the BKY (2006) sharpened q-values;

gen bky06_qval = 1 if pval~=.;

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.;


while `qval' > 0 {;
	* First Stage;
	* Generate the adjusted first stage q level we are testing: q' = q/1+q;
	local qval_adj = `qval'/(1+`qval');
	* Generate value q'*r/M;
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q'*r/M;
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank1 = reject_temp1*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected1 = max(reject_rank1);

	* Second Stage;
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0);
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]));
	* Generate value q_2st*r/M;
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q_2st*r/M;
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank2 = reject_temp2*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected2 = max(reject_rank2);

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition;
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.;
	* Reduce q by 0.001 and repeat loop;
	drop fdr_temp* reject_temp* reject_rank* total_rejected*;
	local qval = `qval' - .001;
};
	

quietly sort original_sorting_order;
pause off;

display "Code has completed.";
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'";
display	"Sorting order is the same as the original vector of p-values";

local outcome "gapman";	
rename bky06_qval `outcome'_qval;
keep `outcome'_qval	;
	tempfile `outcome'_qval;
	save ``outcome'_qval';
	clear;

/*-----------------------------------------------------------------------------------
			 		OUTCOME = Pesticide Residue and Compliance (Table 3)
-----------------------------------------------------------------------------------*/	

	use `maindata';

	local pestoutcome "mean_std_pest_eu comply_china comply_japan comply_eu comply_us";
	foreach p of local pestoutcome {;
	
	preserve;
	keep if round == 2;
	
	local outcome "`p'";
	local count_outcome = `:word count `outcome'';
	mat X = J(7*`count_outcome',1,0);
	local i = 1;
	
	foreach y of local outcome {;
	reg `y' tr_farm tr_int tr_joint certelig tr_fXcert tr_iXcert tr_jXcert `balance' `intbalance' i.strata if round == 2, robust cluster(group_id);

	foreach z of local treat_coef {;
	local t_score = (_b[`z']/_se[`z']);
	local p_score = (2*ttail(e(df_r),abs(`t_score')));
	mat X[`i',1] = `p_score';
	local i = `i' + 1;
	};
	};
	
	clear;
	
*** Anderson's Code for FDR q-value;	
	
	svmat X;
	rename X1 pval;
	
* Collect the total number of p-values tested;

quietly sum pval;
local totalpvals = r(N);

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank;

quietly gen int original_sorting_order = _n;
quietly sort pval;
quietly gen int rank = _n if pval~=.;

* Set the initial counter to 1 ;

local qval = 1;

* Generate the variable that will contain the BKY (2006) sharpened q-values;

gen bky06_qval = 1 if pval~=.;

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.;


while `qval' > 0 {;
	* First Stage;
	* Generate the adjusted first stage q level we are testing: q' = q/1+q;
	local qval_adj = `qval'/(1+`qval');
	* Generate value q'*r/M;
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q'*r/M;
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank1 = reject_temp1*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected1 = max(reject_rank1);

	* Second Stage;
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0);
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]));
	* Generate value q_2st*r/M;
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q_2st*r/M;
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank2 = reject_temp2*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected2 = max(reject_rank2);

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition;
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.;
	* Reduce q by 0.001 and repeat loop;
	drop fdr_temp* reject_temp* reject_rank* total_rejected*;
	local qval = `qval' - .001;
};
	

quietly sort original_sorting_order;
pause off;

display "Code has completed.";
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'";
display	"Sorting order is the same as the original vector of p-values";
	
rename bky06_qval `outcome'_qval;
keep `outcome'_qval	;
	tempfile `outcome'_qval;
	save ``outcome'_qval';
	
	restore;
	};
clear;


/*-----------------------------------------------------------------------------------
			 		OUTCOME = REVENUE AND PROFIT (Table 4, Online Appendix Table A-11)
-----------------------------------------------------------------------------------*/

	use `maindata';
	set seed 1357;
	gen randnum = runiform();
	drop if price == .;
	drop if totalvolume == .;
	drop if cost_costtotal == .;
	drop if revenue == .;
	
	local farmoutcome "log_price log_totalvolume revenue imp_revenue profit imp_profit";
	foreach p of local farmoutcome {;
	
	preserve;
		keep if round == 1;
		sort `p' randnum;
		drop if _n <= 10 | _n > _N-10;
		
	
	
	local outcome "`p'";
	local count_outcome = `:word count `outcome'';
	mat X = J(7*`count_outcome',1,0);
	local i = 1;
	
	foreach y of local outcome {;
	reg `y' tr_farm tr_int tr_joint certelig tr_fXcert tr_iXcert tr_jXcert `balance' `intbalance' i.strata if round==1, robust cluster(group_id);

	foreach z of local treat_coef {;
	local t_score = (_b[`z']/_se[`z']);
	local p_score = (2*ttail(e(df_r),abs(`t_score')));
	mat X[`i',1] = `p_score';
	local i = `i' + 1;
	};
	};
	
	clear;
	
*** Anderson's Code for FDR q-value;	
	
	svmat X;
	rename X1 pval;
	
* Collect the total number of p-values tested;

quietly sum pval;
local totalpvals = r(N);

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank;

quietly gen int original_sorting_order = _n;
quietly sort pval;
quietly gen int rank = _n if pval~=.;

* Set the initial counter to 1 ;

local qval = 1;

* Generate the variable that will contain the BKY (2006) sharpened q-values;

gen bky06_qval = 1 if pval~=.;

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.;


while `qval' > 0 {;
	* First Stage;
	* Generate the adjusted first stage q level we are testing: q' = q/1+q;
	local qval_adj = `qval'/(1+`qval');
	* Generate value q'*r/M;
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q'*r/M;
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank1 = reject_temp1*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected1 = max(reject_rank1);

	* Second Stage;
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0);
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]));
	* Generate value q_2st*r/M;
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q_2st*r/M;
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank2 = reject_temp2*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected2 = max(reject_rank2);

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition;
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.;
	* Reduce q by 0.001 and repeat loop;
	drop fdr_temp* reject_temp* reject_rank* total_rejected*;
	local qval = `qval' - .001;
};
	

quietly sort original_sorting_order;
pause off;

display "Code has completed.";
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'";
display	"Sorting order is the same as the original vector of p-values";
	
rename bky06_qval `outcome'_r1_qval;
keep `outcome'_r1_qval	;
	tempfile `outcome'_r1_qval;
	save ``outcome'_r1_qval';
	restore;	
	
	};		
	
local farmoutcome "log_price log_totalvolume revenue imp_revenue profit imp_profit";
	foreach p of local farmoutcome {;
	
	preserve;
		keep if round == 2;
		sort `p' randnum;
		drop if _n <= 10 | _n > _N-10;	
	
	local outcome "`p'";
	local count_outcome = `:word count `outcome'';
	mat X = J(7*`count_outcome',1,0);
	local i = 1;
	
	foreach y of local outcome {;
	reg `y' tr_farm tr_int tr_joint certelig tr_fXcert tr_iXcert tr_jXcert `balance' `intbalance' i.strata if round==2, robust cluster(group_id);

	foreach z of local treat_coef {;
	local t_score = (_b[`z']/_se[`z']);
	local p_score = (2*ttail(e(df_r),abs(`t_score')));
	mat X[`i',1] = `p_score';
	local i = `i' + 1;
	};
	};
	
	clear;
	
*** Anderson's Code for FDR q-value;	
	
	svmat X;
	rename X1 pval;
	
* Collect the total number of p-values tested;

quietly sum pval;
local totalpvals = r(N);

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank;

quietly gen int original_sorting_order = _n;
quietly sort pval;
quietly gen int rank = _n if pval~=.;

* Set the initial counter to 1 ;

local qval = 1;

* Generate the variable that will contain the BKY (2006) sharpened q-values;

gen bky06_qval = 1 if pval~=.;

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.;


while `qval' > 0 {;
	* First Stage;
	* Generate the adjusted first stage q level we are testing: q' = q/1+q;
	local qval_adj = `qval'/(1+`qval');
	* Generate value q'*r/M;
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q'*r/M;
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank1 = reject_temp1*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected1 = max(reject_rank1);

	* Second Stage;
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0);
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]));
	* Generate value q_2st*r/M;
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q_2st*r/M;
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank2 = reject_temp2*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected2 = max(reject_rank2);

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition;
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.;
	* Reduce q by 0.001 and repeat loop;
	drop fdr_temp* reject_temp* reject_rank* total_rejected*;
	local qval = `qval' - .001;
};
	

quietly sort original_sorting_order;
pause off;

display "Code has completed.";
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'";
display	"Sorting order is the same as the original vector of p-values";
	
rename bky06_qval `outcome'_r2_qval;
keep `outcome'_r2_qval	;
	tempfile `outcome'_r2_qval;
	save ``outcome'_r2_qval';
	restore;	
	
	};		

clear;


/*-----------------------------------------------------------------------------------
			 Trust and Dictator game outcomes (Table 5)
-----------------------------------------------------------------------------------*/

*** Append with Intermediary data;
	use `maindata', clear;
	append using `intdata';
	gen newid = farmerid if farmerid !=.;
	replace newid = traderid+20000 if traderid !=.;
	order newid;
	keep if round == 1;
	gen farmer = (farmerid!=.);


** Merge with game records;
	merge m:1 newid using `game';
	keep if _merge== 3;
	drop _merge;


** Generate variables;

gen g1_partner_int = (g1_partnernewid >= 20000 & g1_partnernewid !=.);
gen g2_partner_int = (g2_partnernewid >= 20000 & g2_partnernewid !=.);
gen g3_partner_int = (g3_partnernewid >= 20000 & g3_partnernewid !=.);	

gen g1_partnershare = g1_partner/100;
gen g2_partnershare = g2_partner/100;
gen g3_partnershare = g3_partner/150;

local g1_partner_char g1_partner_age g1_partner_size g1_partner_volume g1_partner_score g1_partner_type;
local g2_partner_char g2_partner_age g2_partner_size g2_partner_volume g2_partner_score g2_partner_type;
local g3_partner_char g3_partner_age g3_partner_size g3_partner_volume g3_partner_score g3_partner_type;

local treat_coef tr_farm tr_int tr_joint certelig tr_fXcert tr_iXcert tr_jXcert;


	local gameoutcome "g2 g3 g1";
	foreach p of local gameoutcome {;

	preserve;
	keep if farmer == 1;
	
	local outcome "`p'";
	local count_outcome = `:word count `outcome'';
	mat X = J(7*`count_outcome',1,0);
	local i = 1;
	

		foreach y of local outcome {;
	
		qui reg `y'_partnershare `treat_coef' `balance'  i.strata if `y'_partner_int == 1, robust cluster(group_id); 	
		

	foreach z of local treat_coef {;
	local t_score = (_b[`z']/_se[`z']);
	local p_score = (2*ttail(e(df_r),abs(`t_score')));
	mat X[`i',1] = `p_score';
	local i = `i' + 1;
	};
	};
	
	clear;
	
*** Anderson's Code for FDR q-value;	
	
	svmat X;
	rename X1 pval;
	
* Collect the total number of p-values tested;

quietly sum pval;
local totalpvals = r(N);

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank;

quietly gen int original_sorting_order = _n;
quietly sort pval;
quietly gen int rank = _n if pval~=.;

* Set the initial counter to 1 ;

local qval = 1;

* Generate the variable that will contain the BKY (2006) sharpened q-values;

gen bky06_qval = 1 if pval~=.;

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.;


while `qval' > 0 {;
	* First Stage;
	* Generate the adjusted first stage q level we are testing: q' = q/1+q;
	local qval_adj = `qval'/(1+`qval');
	* Generate value q'*r/M;
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q'*r/M;
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank1 = reject_temp1*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected1 = max(reject_rank1);

	* Second Stage;
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0);
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]));
	* Generate value q_2st*r/M;
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q_2st*r/M;
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank2 = reject_temp2*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected2 = max(reject_rank2);

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition;
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.;
	* Reduce q by 0.001 and repeat loop;
	drop fdr_temp* reject_temp* reject_rank* total_rejected*;
	local qval = `qval' - .001;
};
	

quietly sort original_sorting_order;
pause off;

display "Code has completed.";
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'";
display	"Sorting order is the same as the original vector of p-values";
	
rename bky06_qval `outcome'_qval;
keep `outcome'_qval	;
	tempfile `outcome'_qval;
	save ``outcome'_qval';
	restore;

};

clear;

*********** Intermediary game outcomes ************

*** Append with Intermediary data;
	use `maindata', clear;
	append using `intdata';
	gen newid = farmerid if farmerid !=.;
	replace newid = traderid+20000 if traderid !=.;
	order newid;
	keep if round == 1;
	gen farmer = (farmerid!=.);


** Merge with game records;
	merge m:1 newid using `game';
	keep if _merge== 3;
	drop _merge;


** Generate variables;

gen g1_partner_int = (g1_partnernewid >= 20000 & g1_partnernewid !=.);
gen g2_partner_int = (g2_partnernewid >= 20000 & g2_partnernewid !=.);
gen g3_partner_int = (g3_partnernewid >= 20000 & g3_partnernewid !=.);	

gen g1_partnershare = g1_partner/100;
gen g2_partnershare = g2_partner/100;
gen g3_partnershare = g3_partner/150;

local g1_partner_char g1_partner_age g1_partner_size g1_partner_volume g1_partner_score g1_partner_type;
local g2_partner_char g2_partner_age g2_partner_size g2_partner_volume g2_partner_score g2_partner_type;
local g3_partner_char g3_partner_age g3_partner_size g3_partner_volume g3_partner_score g3_partner_type;

local treat_coef tr_farm tr_int tr_joint certelig tr_fXcert tr_iXcert tr_jXcert;

local int_balance firmage dum_firmtype1 dum_firmtype2 dum_firmtype3 size_facility buy_avgprice_fu0 sell_avgprice_fu0 ghp_score_fu0 sell_contract_fu0 sell_exp_fu0 volume_china_fu0 volume_asia_fu0 volume_eu_fu0 volume_dom_fu0;


	local gameoutcome "g2 g3 g1";
	foreach p of local gameoutcome {;

	preserve;
	keep if farmer == 0;
	
	local outcome "`p'";
	local count_outcome = `:word count `outcome'';
	mat X = J(7*`count_outcome',1,0);
	local i = 1;
	

		foreach y of local outcome {;
	
		qui reg `y'_partnershare `treat_coef' `int_balance'  i.strata if `y'_partner_int == 0, robust cluster(group_id); 	
		

	foreach z of local treat_coef {;
	local t_score = (_b[`z']/_se[`z']);
	local p_score = (2*ttail(e(df_r),abs(`t_score')));
	mat X[`i',1] = `p_score';
	local i = `i' + 1;
	};
	};
	
	clear;
	
*** Anderson's Code for FDR q-value;	
	
	svmat X;
	rename X1 pval;
	
* Collect the total number of p-values tested;

quietly sum pval;
local totalpvals = r(N);

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank;

quietly gen int original_sorting_order = _n;
quietly sort pval;
quietly gen int rank = _n if pval~=.;

* Set the initial counter to 1 ;

local qval = 1;

* Generate the variable that will contain the BKY (2006) sharpened q-values;

gen bky06_qval = 1 if pval~=.;

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.;


while `qval' > 0 {;
	* First Stage;
	* Generate the adjusted first stage q level we are testing: q' = q/1+q;
	local qval_adj = `qval'/(1+`qval');
	* Generate value q'*r/M;
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q'*r/M;
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank1 = reject_temp1*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected1 = max(reject_rank1);

	* Second Stage;
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0);
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]));
	* Generate value q_2st*r/M;
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q_2st*r/M;
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank2 = reject_temp2*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected2 = max(reject_rank2);

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition;
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.;
	* Reduce q by 0.001 and repeat loop;
	drop fdr_temp* reject_temp* reject_rank* total_rejected*;
	local qval = `qval' - .001;
};
	

quietly sort original_sorting_order;
pause off;

display "Code has completed.";
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'";
display	"Sorting order is the same as the original vector of p-values";
	
rename bky06_qval `outcome'_int_qval;
keep `outcome'_int_qval	;
	tempfile `outcome'_int_qval;
	save ``outcome'_int_qval';
	restore;

};

clear;


/*-----------------------------------------------------------------------------------
			 		SALES TRANSACTION LEVEL ANALYSIS
-----------------------------------------------------------------------------------*/	
/*-----------------------------------------------------------------------------------
			 		IMPACT ON TRADE WITH INTERMEDIARIES (Table 6)
-----------------------------------------------------------------------------------*/	


*** Read data;
	
	use "${data_path}/FARMER_COMBINETRADE_CLEAN", replace;

	rename followup round;

	merge m:1 farmerid round using "${data_path}/FARMER_LONG_SETUP";
	keep if _merge == 3;
	drop _merge;
	
*** Drop control communes;

	drop if control_commune == 1;	
	
*** Generate dummy for treatment groups;	
	
	tab st_tgroup, gen(dum_tgroup);
	
*** Cluster for standard errors;

	egen group_id = group(strata st_tgroup);

*** Drop large farms (farm size larger than five hectares);

	drop if ha_whitedragon >= 5;	
	
*** Generate treatment levels and interaction terms;

	gen tr_farm = (st_tgroup == 2 | st_tgroup == 6);
	gen tr_int = (st_tgroup == 3 | st_tgroup == 7);
	gen tr_joint = (st_tgroup == 4 | st_tgroup == 8);
	gen certelig = (st_tgroup == 5 | st_tgroup == 6 | st_tgroup == 7 | st_tgroup == 8);
	gen tr_fXcert = tr_farm * certelig;
	gen tr_iXcert = tr_int * certelig;
	gen tr_jXcert = tr_joint * certelig;
	gen iv_tr_farmer = (tr_farm == 1 | tr_joint == 1);
	
** Generate outcome variables;

	gen log_buyer_price = log(buyer_price);
	tab buyer_type, gen(dum_type);
	tab buyer_contract, gen(dum_contract);
	gen nonprogmeet = (1-progmeet);
	gen contract_within = (buyer_contract == 1 | buyer_contract == 2)*progmeet ;
	gen contract_outside = (buyer_contract == 1 | buyer_contract == 2)*nonprogmeet;
	gen spot_within = (buyer_contract == 3)*progmeet ;
	gen spot_outside = (buyer_contract == 3)*nonprogmeet ;
	replace contract_within = 0 if buyer_type == 1 ;
	replace contract_outside = 0 if buyer_type == 1 ;	
	replace spot_within = (buyer_type == 1)*progmeet ;
	replace spot_outside = (buyer_type == 1)*nonprogmeet ;
	
	local balance age female edu_second exp_dragon ha_whitedragon white_age_less1 white_age_less3 white_age_less10 white_age_less20 white_age_more20 anycert_base loan_any_base save_bank_base mean_trust_base mean_business_base mean_confidence_base raven_total present_bias_base gap_comply_base hours_base totalvolume_base price_base log_cost_fert_base log_cost_pest_base log_cost_facequip_base log_cost_labor_base log_cost_utility_base;
	local intbalance int_age int_size_facility int_volume_china_fu0 int_volume_asia_fu0 int_volume_eu_fu0 int_volume_dom_fu0 int_dum_firmtype1 int_dum_firmtype2 int_dum_firmtype3;
	


local contractoutcome "progmeet spot_within spot_outside contract_within contract_outside";
	foreach p of local contractoutcome {;
	
	preserve;
	
	local outcome "`p'";
	local count_outcome = `:word count `outcome'';
	mat X = J(7*`count_outcome',1,0);
	local i = 1;
	
	foreach y of local outcome {;
	qui reg `y' tr_farm tr_int tr_joint certelig tr_fXcert tr_iXcert tr_jXcert `balance' `intbalance' i.strata i.round, robust cluster(group_id);

	foreach z of local treat_coef {;
	local t_score = (_b[`z']/_se[`z']);
	local p_score = (2*ttail(e(df_r),abs(`t_score')));
	mat X[`i',1] = `p_score';
	local i = `i' + 1;
	};
	};
	
	clear;
	
*** Anderson's Code for FDR q-value;	
	
	svmat X;
	rename X1 pval;
	
* Collect the total number of p-values tested;

quietly sum pval;
local totalpvals = r(N);

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank;

quietly gen int original_sorting_order = _n;
quietly sort pval;
quietly gen int rank = _n if pval~=.;

* Set the initial counter to 1 ;

local qval = 1;

* Generate the variable that will contain the BKY (2006) sharpened q-values;

gen bky06_qval = 1 if pval~=.;

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.;


while `qval' > 0 {;
	* First Stage;
	* Generate the adjusted first stage q level we are testing: q' = q/1+q;
	local qval_adj = `qval'/(1+`qval');
	* Generate value q'*r/M;
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q'*r/M;
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank1 = reject_temp1*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected1 = max(reject_rank1);

	* Second Stage;
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0);
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]));
	* Generate value q_2st*r/M;
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals';
	* Generate binary variable checking condition p(r) <= q_2st*r/M;
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.;
	* Generate variable containing p-value ranks for all p-values that meet above condition;
	gen reject_rank2 = reject_temp2*rank;
	* Record the rank of the largest p-value that meets above condition;
	egen total_rejected2 = max(reject_rank2);

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition;
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.;
	* Reduce q by 0.001 and repeat loop;
	drop fdr_temp* reject_temp* reject_rank* total_rejected*;
	local qval = `qval' - .001;
};
	

quietly sort original_sorting_order;
pause off;

display "Code has completed.";
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'";
display	"Sorting order is the same as the original vector of p-values";
	
rename bky06_qval `p'_qval;
keep `p'_qval	;
	tempfile `p'_qval;
	save ``p'_qval';
	restore;	
	
	};	
	
clear;

local contractoutcome "spot_within spot_outside contract_within contract_outside";	
use `progmeet_qval';

	foreach y of local contractoutcome {;
	append using ``y'_qval';
	};	
	
clear;
/*-----------------------------------------------------------------------------------
			 	 APPEND ALL QVAL DATA
-----------------------------------------------------------------------------------*/		


local gapoutcome "gapman";
local pestoutcome "mean_std_pest_eu comply_china comply_japan comply_eu comply_us";
local farmoutcome "log_totalvolume_r1 revenue_r1 imp_revenue_r1 profit_r1 imp_profit_r1 log_price_r2 log_totalvolume_r2 revenue_r2 imp_revenue_r2 profit_r2 imp_profit_r2";
local gameoutcome "g2 g3 g1 g2_int g3_int g1_int";
local contractoutcome "progmeet spot_within spot_outside contract_within contract_outside";	
local appendoutcome "`pestoutcome' `farmoutcome' `gameoutcome' `contractoutcome'";

use `gapman_qval';
foreach y of local appendoutcome {;
	append using ``y'_qval';
	};

save "${data_path}/APPEND_MHT_QVAL.dta"; 


